rm(list=ls())
library(fields); library(maps); library(ncdf4);library(abind)
library(geosphere);library(mapdata)
setwd('/Volumes/Disk2/JJA_ozone_MAM/Harvard_Dataverse')

#----- load functions -----------------------
source("./Function/get_geo.R")
source("./Function/get_met.R")
source('./Function/read_detrend.R')
source('./Function/eof.mca.R')
source('./Function/cov4gappy.R')
source('./Function/Henderson.R')
source('./Function/filled.contour3.R')

#-----------------------------------------
ss1=load("./Figures/Figure 2/mov_detrend.Rdata")

col.gen = colorRampPalette(c("#0563FA", "white","#ff4d4d"), space="rgb")
levels=seq(-0.7, 0.7, length=31)
mid=(levels[1:(length(levels)-1)]+levels[2:(length(levels))])/2
cols=col.gen(length(levels)-1)
if(min(abs(mid))>0){
	ind1=which(mid<=0);ind2=which(mid>=0)
	ind=abind(Reduce(intersect, list(ind1+1, ind2)), Reduce(intersect, list(ind1, ind2-1)))
	ind=sort(unique(ind))
  	 cols[ind]="white"
}	

US.mask=!is.na(apply(yearly_O3,c(1,2),mean,na.rm=TRUE))
US.area=cal.area(lon.frame,lat.frame)
US.area[!US.mask]=NA

#===========================================================
#===========================================================
mai=c(0.3, 0.7, 0.3, 0.1)
mgp=c(1.8, 0.5, 0)
tcl=-0.3
ps=15
par(mai=mai, mgp=mgp, tcl=tcl, ps=ps)
ind1=(lon.frame>=-100 & lon.frame<=-65)
ind2=(lat.frame>=32.5 & lat.frame<=50)

load('./Figures/Figure 2/Henderson.Rdata')
y1=apply(yearly_O3[ind1,ind2,],c(3),mean,na.rm=TRUE)
y2=apply(predict_O3[ind1,ind2,],c(3),mean,na.rm=TRUE) 
# y1=cal.area_mean(yearly_O3,US.area,ind1,ind2)
# y2=cal.area_mean(predict_O3,US.area,ind1,ind2)
R=round(cor((y1),(y2)),digits=2)
plot(1980:2013,(y1),pch=16,col=1,type='o',xlab='',ylab='ozone anomaly (ppbv)',ylim=c(-7,7))
lines(1980:2013,(y2),pch=16,type='o',col="#F74623")
text(paste('r(HF) =',format(R, nsmall=2),sep=''),x=2010,y=5.8,col="#F74623",cex=0.9)

load('./Figures/Figure 2/mov_detrend.Rdata')
y1=apply(yearly_O3[ind1,ind2,],c(3),mean,na.rm=TRUE)
y2=apply(predict_O3[ind1,ind2,],c(3),mean,na.rm=TRUE) 
R=round(cor((y1),(y2)),digits=2)
text(paste('r(MA) =',format(R, nsmall=2),sep=''),x=2010,y=6.8,col="#0060FA",cex=0.9)
lines(1980:2013,(y2),pch=16,type='o',col="#0060FA")

legend(x=1999,y=-3.5,c('obs (HF)','prediction (HF)','prediction (MA)'),col=c(1,"#F74623","#0060FA"),lty=1,pch=16,bty="n",text.col=c(1,"#F74623","#0060FA"))
title(' Mean JJA MDA8 ozone over eastern US',cex.main=0.9,cex.main=0.95,font.main=1)
